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Abstract 

Mean-field coupled lattice maps are used to approximate the physics of driven threshold systems 
I : - .. _ „.. H_ .. . „. . — . ... ... . .e 

^^ . dynamic instability responsible for generating avalanches. Here we present a method of simulating 

Z^ \ specific frictional weakening effects in a mean field slider block model. This provides a means of 

exploring dynamical effects previously inaccessible to discrete time simulations. This formulation 

also results in Abelian avalanches, where rupture propagation is independent of the failure sequence. 

^ \ The resulting event size distribution is shown to be generated by the boundary crossings of a 

/\ ' stochastic process. This is applied to typical models to explain some commonly observed features. 
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I. INTRODUCTION 

In the study of earthquakes, computational slider block models [jj, |2, S 1^ are often 
used to investigate the origin of magnitude-frequency scaling, a ubiquitous feature of global 
seismicity. Slider block models are most commonly implemented as a coupled lattice map, a 
system of continuous variables interacting in discrete time. Similar driven threshold models 
are used to describe a wide range of systems, such as pinned charge density waves |S] , flux 
lattices in type 11 superconductors [u\, and creeping contact lines j7|. All these systems 
exhibit instabilities which form complex spatiotemporal patterns, usually with a regime 
power law events. 



Driven threshold systems with long-range interactions are common in nature, but difficult 



to investigate computationally. However, the 



ning to yield to analytical understanding 



imiting case of mean-field models are begin- 



y,y,yu 



. Mean field model simulations exhibit 



complex event histories and regimes of behavior, including a power law magnitude-frequency 
relation. Despite their simplicity, mean-field models remain sensitive to the choice of up- 
date rules and details of implementation. This is especially evident in models that impose 
a particular form of frictional weakening, where different modes of behavior appear as the 
strength of the weakening is varied J8[ . 

Frictional weakening refers to the drop in cohesive force with velocity or slip, generat- 
ing the dynamical instability which produces an avalanche. It is an essential component 
in the dynamics of rupture, but requires continuous time dynamics to simulate. This is 
prohibitively expensive, severely limiting the scale of the models one can investigate. Cur- 
rent coupled lattice map models necessarily assume the final state does not depend on the 
weakening law, in contradiction to the earliest findings in computational seismology [ij. 

Some coupled lattice map models attempt to simulate frictional weakening with modified 
update rules, generating new species of models with their own behavioral quirks |l^. This 
form of weakening is rigid in design and qualitatively different from the real phenomenon. 
In addition, this method compounds analytical difficulties involving multiple failures of 
individual blocks in a single step J9|. 

Here we present a method of introducing realistic weakening effects in discrete time simu- 
lation. We have developed techniques of simulation and analysis which encompass arbitrary 
weakening laws under identical rules of evolution. These techniques should be applicable to 
a wide range of driven threshold systems, where more precise behavior of the instability may 
be modeled. This unifies the analysis of previously incompatible models, and provides more 
freedom in numerical simulation. This technique also results in Abelian rupture propaga- 
tion, where the size of a simulated earthquake is uniquely determined from initial conditions, 
leading to a rigorous and implementation independent analysis. We then present this anal- 
ysis for mean field slider block models, and discover a theoretical description of observed 
finite-size effects usually interpreted as a nucleation-type phenomenon. 

We first review the basic details of a mean-field slider block model to demonstrate the 
origin of algorithmic dependence. We then present the technique of 'forced weakening', 
which uses a stress excess function, to account for the macroscopic effects of an arbitrary 



microscopic weakening law. This technique is used to understand the relationships between 
current models, and provides a means to explore dynamical effects previously inaccessible 
to discrete time simulation. Finally, we use this technique in model analysis to understand 
the behavior of mean-field models with new accuracy and generality. 

II. THE NEAR MEAN FIELD MODEL 

The general slider-block model represents stick-slip motion along a fault plane with A^ ^ 1 
discrete coordinates (or 'sites') coupled by springs. Each site is assigned a slip deficit Ui which 
measures the distance from global elastic equilibrium. The sites are pinned in place by 
frictional forces, and are subject to a restoring force (which is traditionally called 'stress') 
proportional to their slip deficit. All sites are subject to an external driving force which 
uniformly increases the slip deficits. Internal disorder gives rise to an additional component 
of stress through site interactions. The stress Si at a site i is related to the slip deficits 
through a linear constitutive relation 

Si = -KLUi - ^ Kij {ui - Uj) (1) 

j 

where Kl and Kij are spring constants. If we impose uniform (mean-field) interactions 
between all the elements, Kij = Kc/N, the above relation becomes 

S., = -Klu, - Kciui - {u)) (2) 

where {u) = N~'^ '^i'^i "will denote an average over all A^ sites in the model. We obtain 
a unitless expression by dividing by Kqcl, where a is a characteristic microscopic length 
(the equilibrium distance between sites). Defining the unitless slip deficit = u/a, stress 
a = S/{Kca), and spring constant ratio Kr = K^/Kc, (0) simplifies to 

ai = (0) - (Kn + 1)0,. (3) 

For finite A^ we will refer to this as the near mean-field model. Note that it is easy to invert 
dnj for the slip deficits in terms of stresses, 

, -o-i (a) 



Kr + 1 Kn{KR + 1) 



(4) 



so the configuration is uniquely determined by the parameter Kr and either the shp deficits 
or stresses alone. 

The model is slowly driven away from equilibrium by uniformly increasing all slip deficits. 
Eventually the stress at one site will surpass the maximum local frictional force and 'fail', 
sliding toward its equilibrium point. The motion of a failed site will change the mean slip 
deficit (0), and produce a change in stress at other sites. If this change brings other sites to 
their threshold, they will also fail, producing an avalanche interpreted as a single event. 

This description assumes that when a site begins to slip, the frictional force weakens, 
producing a transient dynamic instability. In discrete time we cannot model the dynamic 
slip or velocity of the site, but instead assign a residual stress a^ at which the motion arrests. 
This 0"^ is chosen from a probability distribution independently for each failed site. Since 
slips occur instantaneously, we lose the interplay between a continuously evolving stress field 
and frictional force at a site. The behavior of dynamical models is known to strongly depend 
on the form of frictional weakening l|, a feature entirely absent from coupled lattice maps. 

Since the near mean field model is not dynamical we are only interested in large-scale 
features of its behavior that are independent of microscopic dynamics. Thus we are free to 
choose the simplest update rules that are consistent with fracture processes. In practice, we 
assume that a single site j reaches its stress threshold first. Since the physics will depend only 
on changes in stress, we may impose a uniform failure threshold a^ by absorbing threshold 
variations into the residual stress distribution. The slip displacement Aj = (jr^ — (jr^ is 
related to the change in stress Actj = a J — a^ by Aj = —Aaj/{Kfi + 1 — N^^). The motion 
of the site will change the mean slip deficit {</)) by Aj/N. We may interpret this as a transfer 
of stress from failing sites to pinned sites. 

A. Model Behavior 

The behavior of this model is observed through numerical simulation and depends on the 
constant Kji, the distribution of residual stresses a^, the initial conditions, and possibly on 
special weakening rules. For large Kji the coupling becomes unimportant and the system 
acts as A^ independent stick-slip blocks. With small Kji and generic (randomized) initial 
conditions the model exhibits a power law in event sizes with the mean-field exponent of 
3/2 (Fig. Q). It is this apparently critical behavior which has drawn attention to this model 
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FIG. 1: Distribution of event sizes for near mean field model simulation with Kji = 0. Note the 
overabundance of large events. 



as an analogue to earthquakes and other largely scale-invariant phenomena. 

In critical models the events are uncorrelated in time or magnitude. Events in any 
magnitude range occur as a Poisson point process in time with a rate appropriate to their 
relative abundance [l6| . The stresses in the model at any time appear uniformly distributed 
between the upper bound of the residual stress and the failure threshold. An example stress 
distribution is shown as a histogram in Fig. |21 

If we assume the stress value of each site is independently sampled from such a uniform 
distribution, and order the results by increasing stress, we expect the stress gaps between 
nearest values to have an exponential distribution. Comparison of this theoretical distri- 
bution with an actual simulation (Fig. ^ show the assumption remarkably valid for large 
N. As long as there is a very small (Order A^~^) randomized residual stress, these stress 
statistics will be persistent in time and independent for each event. With no randomizing 
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FIG. 2: A histogram of stress values of sites in a typical near mean field simulation with uniform 
random residual stress between -0.05 and 0.05. The distribution has the statistical properties of a 
set of uniformly distributed random variables between 0.05 and 1. 



ingredient, limit cycles with time correlated distributions are possible. 

Many ideas concerning the overabundance of large events in Fig. |21have been put forward. 
This has often been considered an example of characteristic behavior, that is, repetition of 
an event whose size is determined by the local geometry. This can be considered a type of 
nucleation phenomena, a runaway Griffith rupture, where events larger than a critical size 
can not arrest. 

This characteristic behavior is emphasized in dynamic weakening models |l5[, where the 
increased propensity for a failed site to slip further (due to a weakened pinning force) is 
simulated by imposing a lower threshold stress cr^ < a^ for the duration of a single event. 
After a site fails it will receive stress transfer from subsequent failures, and thus may reach 
this lower threshold and fail again. Re-failing sites will contribute more to the external 




stress Gap (units of a /N) 

FIG. 3: Comparison of simulated stress gap distribution with expectation from assuming indepen- 
dent uniformly distributed values. Plotted is the cumulative distribution of stress gaps (distance 
between nearest stress values) vs. cumulative exponential distribution. 



transfer and enhance the hkehhood of continued rupture growth. 

This form of weakening manifests itself after some failed sites have their stress brought 
back up to the dynamical threshold. This will only occur when the rupture reaches a certain 
size. After the onset of this dynamical weakening, the additional stress transfer typically 
results in a runaway event encompassing every block in the system. Thus, this form of 
weakening tends to produce characteristic events which always occur once the minimum 
rupture size is reached. To understand how the stress distribution and weakening rules 
determine the event size distribution and model behavior, we must examine the details of 
stress transfer between sites. 



B. Stress Transfer 

The above describes a series type dislocation where sites fail in sequence and the stress 
transfer occurs to other sites all at once. This makes it likely that any failing site (other than 
the single initiator) will have a stress slightly above the threshold, which subtly provides an 
order-of-failure dependence to the stress transfer. As a consequence, the exact stress transfer 
in simulation will depend on obscure factors like the order of iteration over sites. To eliminate 
this we must examine the stress transfer in more detail. 

Suppose that in the course of an event there have been k block failures. Let {k} represent 
the set of indices of failed sites. Call k, = k/N the fraction of failed sites. Then from Q the 
change in stress for any stable site i is 



5_ 

N 
ie{fc} je{k} 



^^^tw = ^ E ^- = ^ E (^i - ^f' 



where 6 = {Kr + 1 — N^^)^^, and {■)j^ is an average applied over failed sites. We call this 
the external stress transfer to signify that it applies to sites that are not part of the rupture. 
The quantity cr| is the stress of site j at failure, which may be greater than a^ . Note that 
the slip displacement Aj is only dependent on the stress drop at failure because pinning 
occurs immediately, and subsequent stress changes will not cause this site to slip further. 
Since a-^ is typically very near a^ and the a^ are identically distributed random variables, 
the term in parenthesis is on average independent of k. Thus the external transfer grows 
linearly with the fraction of failed sites. 

To examine the effects of stress overshoot, suppose a site j slips on the initial step. The 
external transfer is 

(^ij^j = ^i + j;^ (oj 

where for emphasis the superscript denotes the initial stress at a site before any failures 
have occurred. Of course, for the rupture initiator, cr° = cr-^ = a^ . Now suppose this transfer 
causes sites m and n to fail. At this time am,n = o'm,n + ^ {'^j ~ '^f) 1^ i so after m and n 
fail the external transfer becomes 

""^^{^^rnM = ^,, + J^ + J^ + J^ + J^^ . (7) 



Note the factor of 2 in the last term above comes from the fact that two sites {m, n) failed 
in the previous iteration. This simple stress transfer rule introduces dependence on which 
blocks fail during which iteration step. We would prefer the time evolution to be a transpar- 
ent function of the initial stresses only. To satisfy this condition, we must make the stress 
transfer Abelian, that is, independent of the order of failure. This is similar in concept to 
the Abe.,a. .a.d pUe „,ode. QQ. The .ey to ,„aM.g ..ess .ansfe. AbeUan is to use 
the additional degree of freedom offered with the inclusion of simulated weakening. 

C. Forced Weakening 

One way to visualize the effects of weakening is to examine the average stress of sites 
that have failed as a function of rupture size. After failure, a site has average stress (o"'^). 
Subsequent failures will transfer some stress to this now pinned site. Without weakening, 
the average stress of failed sites will grow with rupture size k as 

a-t{K) = ^-{K~N-'){^a),. (8) 

With dynamic weakening, all failed sites with stress > a^ will fail again, putting a ceiling 
on the average internal stress, as illustrated in Fig. 0] 

We might imagine a different approach to weakening where failed sites shed a certain 
fraction of the stress they receive after failure by continuing to slide. This would not be 
convenient to simulate, but clearly possible. This would produce an average internal stress 
function the same as one for larger Kr, however now the stress difference is not dissipated, 
but transferred to external sites. The main point of this would be to generate a weakening 
effect present for all event sizes. Implementing this 'fractional weakening' would involve 
re-computing the slip of all failed sites with each new failure. Instead, we might seek to 
formulate the model so the average internal stress function is given as a physical parameter 
and the requisite slips and stress transfer computed as a result. In essence, in place of 
deriving a weakening function from dynamical equations involving slip or velocity dependent 
friction, we can impose the effects of the weakening as they appear to stable sites. We call 
this approach forced weakening. 

To arrive at a forced weakening formulation we first observe the change in stress of a site 
j as it depends on fc, the number of failed sites. Let A^ denote the slip displacement of 
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FIG. 4: Schematic of average internal stress (average stress of sites that have failed) as a function 
of rupture size for several methods of simulated weakening. 

site j after k failures. If it is nonzero it includes all block motion, including initial failure, 
additional failures from weakening, or continuous sliding. Consider the system of equations 
for the stress changes of the failed sites j G {/c} 

1 v^ 



^^iG{fc} 






(9) 



je{fc} 



where a^ denotes the stress after k failures (and a^ is the initial value). The last line demon- 
strates the simple linear form of the relationship between stress drops and slip displacements. 
This may be obtained via a matrix with diagonal elements A^~^ — {Kr + 1) and off-diagonal 
elements A^~^. This matrix is easily inverted to obtain the slips A^ in terms of the current 
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stress drops Actj 

^^ _ -{Kn + 1 - k/N)Aa^ - N'^ Z^e{k} A^f 
' {KR + l)iKn+l-k/N) 



-M (k/N) {Aa'^)^ ^^^^ 



Kr + 1 {Kr + 1) {Kr + 1- k/N) 

This expression provides the shps that are necessary to generate a set of stress drops. 

The shps are directly related to the external transfer, as in (0). Summing over them 
yields a new expression for the external transfer in terms of the stress drops. 



Ao-,- 



ii{k} 



Ly A, = -(-Aa')^ 



N ^^ ' N 

J6{fc} 



X 



k/N 



Kr + 1 {KR + l){KR+l-k/N)\ 



K 



This is identical to (0) when k = 1. However, the new effective transfer coefficient 6{k,) = 
{Kr + 1 — k)~^ grows with the rupture size. Making up for this is the fact that the stress 
drops are no longer computed only at failure, but account for stress changes occurring as 
the rupture progresses. As failed sites pin and acquire additional stress transfer, the average 
change in stress will decrease. Observe that to calculate the external transfer, we need not 
specify individual stress drops or slips, but require only the average dynamic stress drop of 
all failed sites. This evolving average stress drop is defined as 



(-^-^=lE(-?-- 



k ^'"^ "'- 



= {a')^-a-tiK) (12) 

where we have defined the average internal stress (AIS) function o^{k) = (c^),, which 
characterizes the frictional weakening. The stress transfer now depends only on the average 
initial stress of failed sites (which decreases with rupture size) and the AIS function. Both 
may be defined in terms of rupture size independent of a particular update algorithm. Using 
this formulation, avalanches are Abelian. 

In numerical simulation, we must eventually assign an actual slip and/or residual stress 
to each failed site. Care must be taken to make results agree with the analogous forward 
simulation. For example, given a linear AIS function /(k) = an + P, we could assign the 

12 



k*^ failed site a random residual stress (with mean f3) plus 2a{k — 1)/N. When using AIS 
functions with no forward equivalent, the method of assigning final stresses must be stated 
explicitly. 

The forced weakening method has two main advantages. Practically, it allows the sim- 
ulation of models with arbitrary weakening characteristics, most of which would not be 
obtainable with modified CA rules. Formally, the model is Abelian, so that the event size is 
a unique function of the initial stress configuration. Using this fact we can seek an expression 
which will determine the event size given adequate information of the initial stresses. 

III. MODEL ANALYSIS 

With an Abelian stress transfer, the final event size depends only on the initial stress 
distribution. We could therefore compute the final event size if we concretely characterize 
the stress distribution. For model solutions we find it convenient to mostly work with the 
stress deficits defined as 

S,(t) = a^- a,{t) (13) 

representing the distance a site i is from failure. A rupture originates at sites where S = 
and advances through sites with progressively larger values of S. Consider the cumulative 
distribution -Ps(x) defined as the fraction of sites with stress deficit S < cr^x (scaled so x 
varies from zero to one). Typically, one would choose the microscopic length a such that 
Kcd = cr^ so that a^ = 1. For a specific stress configuration -Ps(x) wiU be a piece- wise 
continuous function with A^ steps. 

A. Solution for Rupture Size 

To express the stress transfer as function of rupture size, the discrete initial stresses aj 
must be expressed in terms of the stress distribution. Beginning with the forced weakening 
form of the external transfer t{k), 

r(«) = (-A^n . (14) 

where (^—Aa'''^ = k~^ '^ji=ik}(^'j ~ '^^)' ^^ notice that the k^^ site to fail in the rupture has 
a stress deficit of -P^ ^(^^)- "^^^ factor A; — 1 refiects the fact that the first site to fail has 
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a stress deficit of zero, so Pj:{0) = N ^, and P^'^{r] < N ^) is undefined. This relationship 
is illustrated graphically in Fig. El Using this fact we may write the average stress drop in 
(HH) as 



cr 



je{k} 

F 



ie{fe} 



7E4K + -?) 



iG{A:} 



0" 



(T — 



K 



Pj, {t] — e)dri — a. 



int\ 



(15) 



where we have identified the average internal stress function (cr^)^ = Wm^k) and introduced 
the notation e = N~'^. Note that we took no limit in writing the integral, but recognized 
the equivalent of the sum and an integral over a step function. 

Using this form of the average stress drop we may write the external transfer in terms of 
the stress distribution as 



Tin, 



Kn + l 



K 



Ka - KaintiK.) - (T / P^ {t] - e)dri 



(16) 



Using (jT^ and the statistical properties oi P^^{ri), we should be able to calculate the size 
of the next event. At any point during a sequence of failures, a rupture will arrest if every 
(initial) stress deficit in the unfailed region is greater than the current external transfer. 
Referring to Figure El we see that the stress deficit of the next site to fail is P^^{k). Using 
this we define a stress excess v which will provide a criteria for the progress of an event in 
terms of the stress distribution: 



v{k) =T{K)-a^Pj.\K) 



(17) 



When V is positive, a rupture will continue to grow, arresting only if/when v{k,) = 0. 
Assuming a rupture has initiated, setting the stress excess to zero yields an integral equation 
for the final rupture size k: 

1 



Kr+I-k 



Kcr - Kaint(H 



a 



j%^\r,-e)dr, 



-a^P^\K) 







(18) 



To solve this equation it would be convenient to consider k a truly continuous variable. The 
stress distribution could vary continuously if we let A^ ^ oo, or by averaging the stress 
deficits over a short time interval. The latter works if we imagine the system driving and 
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FIG. 5: Relationship between inverse cumulative stress deficit function P^; and event size k. = k/N. 
A rupture must initiate at a site where S = 0, so the stress deficit of the k site to fail is given by 
P^iK-N-^). 

stress transfers to be continuous, so all intermediate values of stress are passed through 
during updates. This form of temporal course-graining has been central to previous analysis 
of slider-block models [9|. 

B. Example Solutions 

It is instructive to examine the stress excess function and the event size solutions when 
we can perform the integral over the step function analytically; that is, when the steps are 
of equal size. In this case the stress deficit values are Si = 0, S2 = 1/iV, S3 = 2/N, etc. 
This will also represent the case where N ^ 00, and the empirical distribution of values 
is identical to the probability distribution of a single variable, and the integrals are taken 
over a continuous variable. However, we will not take this limit so factors of e = A^~^ are 
apparent. In this discrete formulation, the integrals become reducible sums. 
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As we observed, in a system exhibiting power law behavior, the stresses are distributed 
as if they were selected from a uniform distribution. Fore simplicity, we can ignore the small 
variation from the uniform distribution due to random residual stresses, since this anomaly 
will only be encountered for ruptures nearly the size of the system. Thus we take Ps(x) = 1, 
P^ix) = f^, and Ps ^(k) = x- 

1. Perfect Weakening 

We will first examine the rather artificial scenario of perfect weakening, where no addi- 
tional stress beyond the residual noise may be supported by failed sites. This is achieved by 
imposing a constant average internal stress ai^(/t) = 0. 

The physical interpretation of this perfect weakening scenario is much like a democratic 
fiber bundle model p^]. Imagine a cable composed of several parallel axial fibers under 
tension, with a gradually increasing load. In this model the individual fibers have some 
randomly distributed strength at which they will break. When a fiber breaks, the load it 
was bearing is distributed equally among all the remaining fibers. 

In the context of a slider block model, failed blocks do not re-pin during a rupture, 
and continue to slide with stress changes from additional failures. The distribution of fiber 
strengths corresponds to the distribution of stress deficits in our system. The factor oil — k 
in the denominator of the stress transfer is equivalent to the rule of sharing the load only 
among unbroken fibers. This earthquake model adds the possibility of stress dissipation 
through the Kfi parameter. When Kr = the models are equivalent. 

The stress excess function for perfect weakening is shown in Fig. IHl Determining the 
rupture size from (fT8|l yields the equation 

k'^ = {2Kr -e)K (19) 

which has solutions for k = or k, = 2Kr — e. For Kn < e the stress excess is always positive 
and the rupture will run away. For slightly larger values of Kn the stress excess function 
starts negative and crosses the axis a short distance from zero, as stress dissipation kills 
the rupture within the first few failures. After crossing zero the function again increases 
without limit. Larger values of Kji do not produce critical event size distributions so aren't 
of interest. With perfect weakening, ruptures of significant size will never arrest. 
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2. No Weakening 

Now consider a model with a linear AIS function, equivalent to typical slider-block models 
with no weakening: 

aint{K > 0) = d- 



2 
The solution to the rupture size equation is now 

K^ = k{2 + 2Kr - E) (20) 

which has solutions at k = and near k = 2 for small Kji. More important, however, is the 
behavior of the stress excess function below k = 1. Again the critical value of K^ is e: for 
Kr < e, the initial stress excess is positive and remains so up to /« = 1, so again any rupture 
that initiates will run away (the model is super-critical). For K^i > e, the stress excess is 
always negative, and no rupture will form (Fig. ^. 

However, when Kji —>■ e the stress excess function vanishes for all k. In this case, we 
have a critically propagating rupture, with the exact stress necessary to tumble each site in 
succession. Also, notice in Fig. IHlhow closely the stress excess remains to zero when Kpt is of 
order e. In this case small fluctuations in the stress distribution could generate a transient 
positive stress excess, resulting in a finite rupture size. 

3. Dynamic Weakening 

Dynamic weakening can be treated with a combination of the above two cases: as long 
as the external transfer is below the dynamic stress threshold, there is no weakening. Once 
external transfer crosses the dynamic threshold, the average internal stress function becomes 
a constant, resulting in runaway stress excess as in the perfect weakening scenario. The point 
at which this crossover happens is easy to calculate by setting t{k.) = cx^ /a^. For W = 
and Kfi = 0, this yields a crossover point at k = cx^ /a^, which is no surprise. Nonzero 
values of W or K^ result in slightly larger values of the crossover point. 

4- Fractional Weakening 

The fractional weakening AIS function is identical to that of no weakening with a modified 
slope. With no weakening the slope is 6/2, which depends on Kr. Fractional weakening 
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FIG. 6: Ideal stress excess functions for models with perfect weakening and no weakening. 

can effectively balance the effect of dissipation, moving the critical value of K^ up. Other 
conclusions of the case with no weakening remain the same. 

Using these ideal solutions we can easily understand the role of the control parameters and 
weakening function in terms of how a single event propagates. Clearly, without fluctuations 
in the stress distribution, the model behavior is trivial. Now we will examine solutions for 
finite models where the time-averaged stress distribution is a stochastic process. 

IV. FLUCTUATIONS AND STOCHASTIC RUPTURE PROPAGATION 

When using a smooth stress distribution the solutions for the rupture size are trivial, 
producing no complexity. Passing to the thermodynamic limit yields a model with no fluc- 
tuations in the stress distribution, a true mean field model. This is not what we observe in 
simulation, as actual stress distributions are finite realizations. A model with a finite number 
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of elements, no matter how large, is inherently different from an infinite model. The finite 
A^ and consequent course grained empirical distribution provides the inherent discreteness 
necessary to produce complex behavior, a near mean field model. 

In a short time interval At where no failures occur, the stress at each site changes an 
amount Ao" = {Kfi+l)VAt where V is a velocity measured in units of the microscopic length 
a per unit time. A time average is equivalent to specifying -Ps(x) as the integral of a density 
function Pe(x)? which is the average number of sites with stress deficit S between a^x ^^^ 
a^{x + Ax), where A^ is the change in scaled stress (Acr/a^) over the time interval. For a 
stress distribution of normalized width / (0 < / < 1), the expected time to the next event 
is defined by Ax = I/N, or At = I/[NV{Kii + 1)]. 

Based on the mean field distribution of event sizes the expected event size at any time is 
k = yN . Thus, during the average time interval between events, we expect the distribution 
function to average over v^ different stress values. For large N , this should be adequate 
to define a stress distribution function as a stochastic process with well defined mean and 
variance at each point. If we use this statistical description of the stress distribution, we 
will obtain a statistical answer for the event size. 

The near mean field model preserves the essential fluctuations characteristic of discrete 
models, while also proving analytically tractable. The stress distribution averaged over the 
inter-event time is a stochastic process, and computing the event size from this will yield a 
probability distribution. Since we observe the stress distribution is statistically stationary, 
the distribution of event sizes for any time step is also the time-average distribution. 

Again consider the case where the stresses (or stress deficits) appear to have been chosen 
from a uniform distribution. Consider an ordering of stress deficits Si < Sj < . . . < S^r 
and define Si = and S^v+i = 1. Let Yi = Sj+i — Sj for i = 1 . . .N. Then we would 
expect the gaps Y^ to be independent identically distributed (i.i.d.) random variables with 
an exponential distribution of mean fiy = N~'^ and variance Vary = A^~^ (taking W = 
for simplicity). Derivations of the event size distribution from these statistics have been 
considered before [sl, luj, but the forced weakening formulation provides new rigor and 
generality in the results. We will also derive the finite-size correction to the mean field 
power law behavior for the first time. 

Notice that the partial sum over intervals is nothing but the inverse cumulative distribu- 
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tion, 



E^^ 


- Sfc+1 — El 




- ^=-'(^' 



(21) 

iV 

yielding the stress deficit of tlie /c*^ site to fail. 
Next define the zero mean random variables 

= v^Vf.-^ (22) 

that are partial sums of the gaps up to the k^^ site. As A^ ^ oo the intervals Xfi — Xt are 
Gaussian random variables with mean zero and variance t' — t. This is a standard method of 
constructing a Wiener process from any underlying i.i.d. sequence of random variables. Since 
we aren't actually taking the infinite A^ limit, we should consider this process equivalent to 
a Wiener process at scales much greater than A^~^. 

Knowing that the Xt represents a Brownian motion in the continuum limit we may express 
the inverse cumulative distribution as 

.,H.,^E.^(^4 

^'^ +K (23) 



'N 

Thus the inverse cumulative distribution may be presented as a sum of the ideal linear term 
(drift) and a Wiener process with variance A^^^. Note that the fluctuations scale with the 
system size as expected. 

The stress excess function also depends on the integral of the inverse cumulative distri- 
bution, representing the average stress deficit of failed sites 

/.(S), = £p^\x)dx = Y + ^£x^dx (24) 

= Y + ^. (25) 

This integral introduces a new process S,^. This is a sum of independent Gaussian random 
variables with zero mean. The result is a random variable with zero mean and a variance 
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K^/3N + 0{1/N'^) jl6|. This term does not bias the stochastic component, and for large A^ 
it's contribution will be negligible next to X^/yN . 
The stress excess function now becomes 

„2 



v{k) 1 



a^ Kn + l 






K-^. (26) 



N 



First examine the behavior of this stochastic function near k = 0. Consider the critical 
model {Kn = N~^) with no weakening. In studying (J2(J|) we saw this function vanishes 
for all values of k. Only the noise terms, dominated by X^/V^, are present. Thus, the 
distribution of event sizes is the distribution of zero-crossings of a standard random walk. 

A. Event Size Distributions 

The zero-crossing distribution is easily derived using the fact that any i.i.d. distribution 
can be used to construct the Wiener process with identical results. It's convenient to build 
the walk from a discrete random variable taking values of ±1 with equal probability. Let 
Pd{k) denote the probability that the walk is d steps from the origin after k steps are taken. 
Taking the (horizontal) step size of the walk to be 1/2N, returns to zero are possible at 
steps k/N with probability 

p(k)- ^ (^^-^)' r27i 

^o(fc) - ^(^_i),2 (27) 

The probability of k being the first crossing of zero is Po{k)/k. These probabilities have a 
power law distribution as seen by taking the log and applying Stirling's approximation: 

pk L-3/2 

^ ^ ^ (28) 

The result is a power law with an exponent of —3/2 as observed in Fig. ^ This distribu- 
tion matches the power law behavior observed in simulation, except at the largest values. 

The expected distance d from the origin for an unbiased random walk in terms of the 
number of steps k is given by {d^) — ^Jhnj2 for large k. With our distance per step of 
1/ \/N , the expected magnitude of the walk grows as ^Jt:k/2. However, all the stress gaps 
must add to one (minus the 0{1/N) distance a^ — Sjv), so the walk must be constrained 
to return to zero at k = 1. If the stress values are considered a Poisson point process, this 
is equivalent to requiring there to be exactly N values between and a^ . For N ^ 1, this 
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constraint is unnoticeable in the statistics of individual stress gaps, but manifests itself when 

all the gaps are summed. 

To correct this we weight each probability Po{k) by the probability that a walk starting 

at that point will return to zero at A; = A^. This latter quantity is Po{N — k). Dividing 

the product of these terms by Po{N) provides the proper normalization. The resulting 

constrained distribution P[^{k) is given by 

_ (2fc)!(2iV-2fe)!iVP 

'^'''~ i2N)\iN-k)m^ ■ ^^^ 

The approximated first passage distribution P^ik) is shown in Fig. [7| This is well approxi- 
mated by a corrected power law 

m^) = r^ n (30) 

Note this has a minimum at exactly /t = 3/4. We conclude that this feature, which resembles 
a nucleation phenomenon, stems from the basic statistics of the stress distribution. For 
models with a large number of sites, the low number of large events often makes it difficult 
to compare the data to an exact power law. However, we can also compare the cumulative 
stress distribution the the integral of (j3Up . This is shown in Fig. |H1 However, a discrepancy 
at small event sizes, due to error in Stirling's approximation is evident in this graph. Using 
the exact combinatorial expression (j^^ instead yields a perfect fit to the simulation data. 

Note in this particular case we can also derive the distribution from first principles, since 
each site contributes exactly 1/A^ to the stress transfer, when the rupture arrests there must 
be exactly k sites in an interval k/N, which for N uniformly distributed variables can be 
written 

However, the above treatment is more general, and may be used to consider the distribution 
of event sizes for non-critical models. 

When Kr > e, The idealized stress excess is negative (and diverges at k = K^ + 1). The 
magnitude of the stochastic term must overcome this deficit for a rupture to propagate. A 
previous approach J8[ is to linearize the stress excess about k = (not including terms of 
0{e)) 

v(k) ^ K ^= (32) 

^ ' 1 + Kr ^ ^ ^ 
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FIG. 7: Theoretical event size distribution for a no-weakening critical model with the additional 
constraint on the stochastic component of the stress distribution. This results in an exact fit to 
the simulation data. 



and compute the first crossings of the random walk with the line. The exact result is given 
by the Bachelier-Levy formula [17 1. We can derive the result by defining a biased walk where 
an step up has probability p = \ — f^ and a step down q = \ + -f^, where c = 2{k\i) - 
The first crossing probabilities for A; = 1, 2, . . . are 



Po{k) l{2k)\ [I c 



^-3/2g-4c2. iV » 1. 



(33) 



k k kP \A N 

For the limiting value of c = 1/2 {Kn — ^ oo), the power law is dominated by an exponential 
decay e"''. 

However, the linearized stress excess does not provide a valid approximate solution. The 
actual distribution will be truncated much more quickly. Ignoring terms of 0{e) we can 
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FIG. 8: The Cumulative distribution of events compared to the theoretical power law. Discrepan- 
cies for small events show the limits of Stirling's approximation. Using the combinatorial expression 
yields an exact fit to simulation data. 



write the entire expression for the stress excess 



iy{K) 



K 



Kr + 1 



-K 



R 



K 



^R 
R 



2{Kr + 1){Kr 



-K 



K 



+ 



K 



Kr+1 ■ Kr + 1~k\ ^^'^^ 

Note that the magnitude of the nonUnear term overtakes the first when k > {Kr + l)/2. 
The second hne gives the stress excess function in its most compact form. 

If we label the non-stochastic part of the stress excess z/(/t) = ^'{k) + Xf^/y/N, the event 
size distribution can be formulated as the distribution of first intersection (Fig. E)). 



Xk 



N \u'{k)\ 



(35) 



where we recognize the sign of the walk term is irrelevant. There is an extensive literature 
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FIG. 9: Event size distributions for simulations with Kji > e. Also shown is the theoretical 
distribution based on the linearized stress excess for Kr = 1. The data falls off much more rapidly 
due to the nonlinear term. 

regarding the intersection of a random walk with a curvilinear boundary [l7|, |l8|. Several 
formal solutions are available, but analytic forms exist only for special cases. Numeric 
computation based on the formal solution is still useful because it eliminates sampling errors 
for low probability events. 

V. CONCLUSIONS 

The forced weakening method introduces a means of modeling dynamical behavior leading 
to an instability in an efficient discrete time simulation. This is done by using an arbitrary 
weakening function to compute stress transfers during simulation, then assigning random 
residual stresses consistent with the weakening function at the completion of rupture. This 
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technique results in Abelian avalanches allowing a rigorous and implementation-independent 
analysis. Event distributions derived from this analysis yield excellent agreement between 
theory and simulation, including a novel explanation for the overabundance of large events. 
If one used a small-scale dynamical model to characterize the average internal stress of a 
rupture as it depends on the physical parameters of a rate and state dependent frictional 
law, we could understand features of the resulting event size distribution. 

This approach should be equally applicable to related discrete threshold models. While 
the resulting formalism is dependent on the mean-field character of the model, the numerical 
techniques may find wider use whenever an inversion like (jlOj) is available. Use of a stress 
excess function to understand the event size distribution may be extendible to non mean 
field models with a finite interaction range. Local correlations might result in the stochas- 
tic component resembling a fractional Brownian motion, predictably altering the crossing 
distribution. 

E.F.P. and J.B.R. acknowledge support by DOE grant DE-FG03-95ER14499; J.S.S.M. 
was supported as a Visiting Fellow by CIRES, University of Colorado at Boulder. 
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